
get_coefs <- function(type, n_covs, R2_d, R2_y) {
  set.seed(02140)
  sigma <- toeplitz(0.5 ^ (0:(n_covs - 1)))

  if (type == "main") {
    b <- (1 / (1:n_covs)) ^ 2
    beta_y <- c(2 * b)
    beta_v <- c(2 * b)
    beta_d <- c(2 * b)

    ## get approx cov matrix
    X <- MASS::mvrnorm(n = 10000, mu = rep(0, n_covs), Sigma = sigma)
    X <- scale(X, scale = FALSE)
    c_names <- paste("X", 1:n_covs, sep = "")
    colnames(X) <- c_names
    # now generate V (moderator) and D (treatment) using X and relevant betas
    V <- rbinom(10000, 1, boot::inv.logit(X %*% beta_v))
    XV <- X * V
    xv_var1 <- var(XV - fitted(lm(XV ~ X + V)))
    kk1 <- c(t(b) %*% xv_var1 %*% b)
    c_d <- c(sqrt(R2_d / ((1 - R2_d) * kk1)))

    beta_vd <- c(c_d * b)
    d <- X %*% beta_d + V * 0.25 + XV %*% beta_vd + rnorm(10000)
    xv_var2 <- var(XV - fitted(lm(XV ~ d * V + X + V)))
    kk2 <- c(t(b) %*% xv_var2 %*% b)
    c_y <- c(sqrt(R2_y / ((1 - R2_y) * kk2)))
    beta_vy <- c_y * b
    beta_dy <- 0 * b
    ##   dX <- c(d) * X
    ##  colnames(dX) <- paste("dX", 1:n_covs, sep = "")
    ## e_y <- X %*% beta_y + d * treatment_effect + V * moderator_effect +
    ##   alpha0 * d * V + (XV) %*% beta_vy + dX %*% beta_dy
    ## y_sd <- 1 ##sqrt((1 + e_y) ^ 2 / mean((1 + e_y)^2) )
    ## y <- e_y  + rnorm(10000, sd = y_sd)

  } else if (type == "binary") {

    b_y <- c(1/c(1:5), rep(0, times = 5), 1/c(1:5), rep(0, times = n_covs - 15))
    b_d <- c(1/c(1:10), rep(0, times = n_covs - 10))
    beta_y <- c(0.1 * b_y)
    beta_v <- c(0.1 * b_d)
    beta_d <- c(0.5 * b_d)

    ## get approx cov matrix
    X <- MASS::mvrnorm(n = 10000, mu = rep(0, n_covs), Sigma = sigma)
    X <- scale(X, scale = FALSE)
    c_names <- paste("X", 1:n_covs, sep = "")
    colnames(X) <- c_names
    # now generate V (moderator) and D (treatment) using X and relevant betas
    V <- rbinom(10000, 1, boot::inv.logit(X %*% beta_v))
    XV <- X * V
    xv_var1 <- var(XV - fitted(lm(XV ~ X + V)))
    kk1 <- c(t(b_d) %*% xv_var1 %*% b_d)
    c_d <- c(sqrt(R2_d / ((1 - R2_d) * kk1)))


    beta_vd <- c(c_d * b_d)
    d <- X %*% beta_d + V * 0.25 + XV %*% beta_vd + rnorm(10000)
    xv_var2 <- var(XV - fitted(lm(XV ~ d * V + X + V)))
    kk2 <- c(t(b_y) %*% xv_var2 %*% b_y)
    c_y <- c(sqrt((pi^2 / 3) * R2_y / ((1 - R2_y) * kk2)))
    beta_vy <- c_y * b_y
    beta_dy <- 0 * b_y

  } else if (type == "dense") {
    
    b <- (1 / (1:n_covs))
    beta_y <- c(2 * b)
    beta_v <- c(2 * b)
    beta_d <- c(2 * b)

    ## get approx cov matrix
    X <- MASS::mvrnorm(n = 10000, mu = rep(0, n_covs), Sigma = sigma)
    X <- scale(X, scale = FALSE)
    c_names <- paste("X", 1:n_covs, sep = "")
    colnames(X) <- c_names
    # now generate V (moderator) and D (treatment) using X and relevant betas
    V <- rbinom(10000, 1, boot::inv.logit(X %*% beta_v))
    XV <- X * V
    xv_var1 <- var(XV - fitted(lm(XV ~ X + V)))
    kk1 <- c(t(b) %*% xv_var1 %*% b)
    c_d <- c(sqrt(R2_d / ((1 - R2_d) * kk1)))

    beta_vd <- c(c_d * b)
    d <- X %*% beta_d + V * 0.25 + XV %*% beta_vd + rnorm(10000)
    xv_var2 <- var(XV - fitted(lm(XV ~ d * V + X + V)))
    kk2 <- c(t(b) %*% xv_var2 %*% b)
    c_y <- c(sqrt(R2_y / ((1 - R2_y) * kk2)))
    beta_vy <- c_y * b
    beta_dy <- 0 * b
    
  } else if (type == "sparse") {
    
    b <- (1 / (1:n_covs)) ^ 2
    b[-c(1:15)] <- 0
    beta_y <- c(2 * b)
    beta_v <- c(2 * b)
    beta_d <- c(2 * b)

    ## get approx cov matrix
    X <- MASS::mvrnorm(n = 10000, mu = rep(0, n_covs), Sigma = sigma)
    X <- scale(X, scale = FALSE)
    c_names <- paste("X", 1:n_covs, sep = "")
    colnames(X) <- c_names
    # now generate V (moderator) and D (treatment) using X and relevant betas
    V <- rbinom(10000, 1, boot::inv.logit(X %*% beta_v))
    XV <- X * V
    xv_var1 <- var(XV - fitted(lm(XV ~ X + V)))
    kk1 <- c(t(b) %*% xv_var1 %*% b)
    c_d <- c(sqrt(R2_d / ((1 - R2_d) * kk1)))

    beta_vd <- c(c_d * b)
    d <- X %*% beta_d + V * 0.25 + XV %*% beta_vd + rnorm(10000)
    xv_var2 <- var(XV - fitted(lm(XV ~ d * V + X + V)))
    kk2 <- c(t(b) %*% xv_var2 %*% b)
    c_y <- c(sqrt(R2_y / ((1 - R2_y) * kk2)))
    beta_vy <- c_y * b
    beta_dy <- 0 * b
    
  }
  return(list(beta_d = beta_d, beta_y = beta_y, beta_v = beta_v,
              beta_vy = beta_vy, beta_vd = beta_vd, beta_dy = beta_dy))
}
